// Copyright 2012 The Chroma Authors. All rights reserved. See the LICENSE file.

// Spectral Power Distribution of a CIE D-Illuminant
// The relative spectral power distribution of a CIE D-Illuminant is computed from its chromaticity coordinates. 
// The chromaticity coordinates may be computed from the correlated color temperature.
// Ref.: [28] 

package white

import "math"

func s(λ float64) (s0, s1, s2 float64) {
	s := [54][3]float64{
		{0.04, 0.02, 0.0},
		{6.00, 4.50, 2.0},
		{29.60, 22.40, 4.0},
		{55.30, 42.00, 8.5},
		{57.30, 40.60, 7.8},
		{61.80, 41.60, 6.7},
		{61.50, 38.00, 5.3},
		{68.80, 43.40, 6.1},
		{63.40, 38.50, 3.0},
		{65.80, 35.00, 1.2},
		{94.80, 43.40, -1.1},
		{104.80, 46.30, -0.5},
		{105.90, 43.90, -0.7},
		{96.80, 37.10, -1.2},
		{113.90, 36.70, -2.6},
		{125.60, 35.90, -2.9},
		{125.50, 32.60, -2.8},
		{121.30, 27.90, -2.6},
		{121.30, 24.30, -2.6},
		{113.50, 20.10, -1.8},
		{113.10, 16.20, -1.5},
		{110.80, 13.20, -1.3},
		{106.50, 8.60, -1.2},
		{108.80, 6.10, -1.0},
		{105.30, 4.20, -0.5},
		{104.40, 1.90, -0.3},
		{100.00, 0.00, 0.0},
		{96.00, -1.60, 0.2},
		{95.10, -3.50, 0.5},
		{89.10, -3.50, 2.1},
		{90.50, -5.80, 3.2},
		{90.30, -7.20, 4.1},
		{88.40, -8.60, 4.7},
		{84.00, -9.50, 5.1},
		{85.10, -10.90, 6.7},
		{81.90, -10.70, 7.3},
		{82.60, -12.00, 8.6},
		{84.90, -14.00, 9.8},
		{81.30, -13.60, 10.2},
		{71.90, -12.00, 8.3},
		{74.30, -13.30, 9.6},
		{76.40, -12.90, 8.5},
		{63.30, -10.60, 7.0},
		{71.70, -11.60, 7.6},
		{77.00, -12.20, 8.0},
		{65.20, -10.20, 6.7},
		{47.70, -7.80, 5.2},
		{68.60, -11.20, 7.4},
		{65.00, -10.40, 6.8},
		{66.00, -10.60, 7.0},
		{61.00, -9.70, 6.4},
		{53.30, -8.30, 5.5},
		{58.90, -9.30, 6.1},
		{61.90, -9.80, 6.5}}

	if λ < 300 || λ > 830 {
		panic("λ out of range")
	}
	i := int(math.Floor((λ - 300) / 10.0))
	s0, s1, s2 = s[i][0], s[i][1], s[i][2]
	return
}

func SD(λ, xD, yD float64) (sd float64) {

	m := 0.0241 + 0.2562*xD - 0.7341*yD
	m1 := (-1.3515 - 1.7703*xD + 5.9114*yD) / m
	m2 := (0.0300 - 31.4424*xD + 30.0717*yD) / m

	t1, t2, t3 := s(λ)
	t2 *= m1
	t3 *= m2
	sd = t1 + t2 + t3
	return
}
